    Info<< "Reading field p\n" << endl;
    volScalarField p
    (
        IOobject
        (
            "p",
            runTime.timeName(),
            mesh,
            IOobject::MUST_READ,
            IOobject::NO_WRITE
        ),
     mesh
    );

    scalar avgRho = 1;

    // determine an average density for the p/T conditions
    {
        autoPtr<basicPsiThermo> thermo
        (
            basicPsiThermo::New(mesh)
        );

        volScalarField rho
        (
            IOobject
            (
                "rho",
                runTime.timeName(),
                mesh,
                IOobject::READ_IF_PRESENT,
                IOobject::AUTO_WRITE
            ),
            thermo->rho()
        );

        avgRho = gAverage(rho);

	Info<< "average density = " << avgRho << nl << endl;
    }

    p = dimensionedScalar("zero", p.dimensions(), 0.0);

    Info<< "Reading field U\n" << endl;
    volVectorField U
    (
        IOobject
        (
            "U",
            runTime.timeName(),
            mesh,
            IOobject::MUST_READ,
            IOobject::AUTO_WRITE
        ),
        mesh
    );

    U = dimensionedVector("0", U.dimensions(), vector::zero);

    // avoid writing phi (dimensional mis-match with other problems)
    surfaceScalarField phi
    (
        IOobject
        (
            "phi",
            runTime.timeName(),
            mesh,
            IOobject::NO_READ,
            IOobject::NO_WRITE
        ),
        fvc::interpolate(U) & mesh.Sf()
    );

    // now that phi has been created, we can impose the velocity bcs
    // we need this since the momentum equations are not solved
    U.correctBoundaryConditions();

    label pRefCell = 0;
    scalar pRefValue = 0.0;
    setRefCell(p, mesh.solutionDict().subDict("SIMPLE"), pRefCell, pRefValue);
